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Three-dimensional topological insulators harbour metallic surface states with exotic properties. 

In transport or optics, these properties are typically masked by defect-induced bulk carriers. Com¬ 
pensation of donors and acceptors reduces the carrier density, but the bulk resistivity remains 
disappointingly small. We show that measurements of the optical conductivity in BiSbTeSe 2 pin¬ 
point the presence of electron-hole puddles in the bulk at low temperatures, which is essential for 
understanding DC bulk transport. The puddles arise from large fluctuations of the Coulomb po¬ 
tential of donors and acceptors, even in the case of full compensation. Surprisingly, the number of 
carriers appearing within puddles drops rapidly with increasing temperature and almost vanishes 
around 40 K. Monte Carlo simulations show that a highly non-linear screening effect arising from 
thermally activated carriers destroys the puddles at a temperature scale set by the Coulomb inter¬ 
action between neighbouring dopants, explaining the experimental observation semi-quantitatively. 
This mechanism remains valid if donors and acceptors do not compensate perfectly. 

PACS numbers: 72.20.-i, 74.62.Dh, 78.20.-e, 78.30.-j 


I. INTRODUCTION 


Three-dimensional topological insulators attract sig¬ 
nificant attention mostly because they feature two- 
dimensional Dirac fermions on the surface that possess 
the peculiar characteristics of spin-momentum locking 
and topological protection m- The existence of such 
Dirac fermions has been confirmed by surface-sensitive 
techniques such as angle-resolved photoelectron spec¬ 
troscopy or scanning-tunneling microscopy [SHin]. How¬ 
ever, the exotic phenomena expected in the electromag¬ 
netic response of these systems largely remain unex¬ 
plored to date. Prominent examples are a topological 
magnetoelectric effect related to the quantum Hall ef¬ 
fect of the surface states yielding magnetic monopoles 
as mirror charges of electric charges mm , a univer¬ 
sal Faraday rotation angle given by the vacuum fine- 
structure constant a diia, or a universal surface con¬ 
ductance G{uj) = TTB^/Sh for energies Huj larger than twice 
the Fermi energy Ep min]. These effects are blurred 
by a dominant bulk conductivity in real specimen of 
topological insulators such as the prototypical binary 
tetradymites Bi 2 Te 3 and Bi 2 Se 3 , which can be catego¬ 
rized as degenerate semiconductors. Typically, single 
crystals of these compounds show defect-induced charge 
carriers with densities above a few 10^® cm’ -3 dMi]. 

Understanding the defect chemistry allowed for a dra¬ 
matic reduction of the carrier density H ng. Near- 
stoichiometric Bi 2 Se 3 exhibits n-type conductivity orig¬ 
inating from Se vacancies acting as donors, whereas p- 
type conductivity predominates in Bi 2 Te 3 due to antisite 
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defects. The most successful route to reduced bulk con¬ 
ductivity aims at two goals in parallel: reduction of the 
defect density and compensation of the remaining de¬ 
fects, i.e., K = Na/Nd = 1, where Np and Na denote 
the densities of donors and acceptors, respectively. In 
Bi 2 _a;Sb 2 ,Te 3 _ySey, a reduced defect density is ac hieved 
by chalcogen order [23l[24] (see Methods, Sec. VIA) while 
variation of x and y allows for optimized compensation 
in combination with the possibility to tune the energy 
of the Dirac point with respect to the Fermi energy Ep 
[25l[26]. In BiSbTeSe 2 , the Dirac point nearly coincides 
with Ep, it thus may serve as a benchmark for the bulk 
carrier dynamics at very low carrier concentrations. 


For a sample thickness d < 10 pm, the bulk conduc¬ 
tance of BiSbTeSe 2 is low enough at low temperatures to 
be out-weighted by the surface conductance [231 [271128]. 
This allows to observe a hallmark of topological trans¬ 
port, the half-integer quantum Hall effect, at tempera¬ 
tures up to 35K [57]. The bulk resistivity pb{T) never¬ 
theless raises several questions. At low temperatures, 
Pb{T) of thick samples of Bi 2 _a;Sba;Te 3 _j,Sey and also of 
Bi 2 Te 2 Se does not exceed 10-20Hem [T71 [75H55] . 

even when the shunting effect of the surface is taken 
into account [28]. Above about 100 K, pb[T) shows acti¬ 
vated behavior !xexp{EA/kBT) but the activation energy 
Ea appears to be substantially smaller than the intrin¬ 
sic value given by half the gap size, A/2 [ 21 ]. The bulk 
conduction mechanism of this important class of materi¬ 
als should be better understood and controlled for future 
investigations of novel topological phenomena. 


A theoretical explanation of the small activation en¬ 
ergy has recently been suggested by Skinner, Chen, 
and Shklovskii [36U38] building on previous work [39]. 
They considered a perfectly compensated semiconduc¬ 
tor (Np = Na = Ndef) with shallow donor and acceptor 
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levels. In such a system, donors give electrons to ac¬ 
ceptors, resulting in positively charged donors and nega¬ 
tively charged acceptors. In this situation, however, the 
long-range Coulomb interactions necessarily enforce the 
formation of large puddles, i.e., regions in the bulk which 
are either p- or n-doped. The reason is that in a vol¬ 
ume of size i?^, random fluctuations of the donor and 
acceptor densities Nd and Na lead to a typical charge 
of order ey/N^efR^ and therefore to a Coulomb poten¬ 
tial of order e'^y/N^efR^/, where £ denotes the 
dielectric constant and e the elementary charge. The po¬ 
tential fluctuations grow proportional to VR and become 
as large as A/2 at a length scale Rg = (A/Ec)'^ ddef/^ir 
[36] which is much larger than the average defect distance 

_T /q 

ddef = -^def^ experimentally relevant case A^ E^, 

where i?c = e^/(47r£o£ddef) denotes the Coulomb interac¬ 
tion between neighbouring dopants. On this length scale 
Rg, the valence and conduction bands are deformed so 
strongly that they touch and cross the chemical potential, 
giving rise to electrically conducting puddles, see Fig. 

In the case of electron puddles, for example, some of the 
donor states become occupied resulting in neutral donors 
in this region. Based on this puddle scenario, Shklovskii 
and coworkers [36] find activated behavior of the resis¬ 
tivity, pb{T) oc exp{E a/ ksT), above roughly 40 K with 
a small activation energy Ea ~ 0.15 A, consistent with 
experimental values |24| . At lower temperatures one ex¬ 
pects to observe the famous Efros-Shklovskii law [3D] for 
variable-range hopping, pb{T) oc exp[(rEs/ir)^/^]. This, 
however, does not describe the experimentally observed 
small bulk resistivity at the lowest temperatures. Never¬ 
theless, the physics of puddle formation is a prime can¬ 
didate to explain why it is so difficult to reach high bulk 
resistivities in compensated topological insulators. 

A direct experimental detection of electrically conduct¬ 
ing puddles in the bulk of topological insulators is there¬ 
fore highly desirable. Surface-sensitive techniques are not 
ideally suited, as puddles are strongly suppressed close 
to the metallic surface which provides an extra screen¬ 
ing channel m- Nevertheless, the size of potential fluc¬ 
tuations observed in scanning tunneling microscopy m 
appears to be consistent with puddle formation [37] . 

Optical spectroscopy is a bulk-sensitive method ide¬ 
ally suited to detect large conducting regions. The op¬ 
tical properties of Bi 2 Te 3 , Bi 2 Se 3 , and of solid solutions 
thereof were investigated intensively already half a cen¬ 
tury ago [ITVHH] in view of their favorable thermoelectric 
properties [46] |47| . Recently, optical data were reported 
for single crystals of Bi 2 Te 2 Se and Bi 2 _ 2 ,Sb 3 ;Te 3 _ySey 
showing reduced carrier density [35H55]. However, the 
bulk carrier dynamics at very low densities were not ad¬ 
dressed in detail. In particular, these data do not allow 
to draw conclusions on the presence of puddles. Here, 
we give a detailed account of the optical properties of 
the approximately fully compensated topological insula¬ 
tor BiSbTeSe 2 in the infrared range. We reveal clear 
signatures of conducting puddles, making use of the re¬ 
cent achievement [24] of very low carrier densities in 



FIG. 1: Illustration of puddle formation. The left panel de¬ 
picts the spatial variation of the energies E±{r) = F(r)±A/2 
of conduction and valence bands (upper and lower lines) 
caused by the long-ranged Coulomb potential V (r) arising 
from randomly placed donors and acceptors. At T = 0, the 
bands fluctuate so strongly that the chemical potential p is 
crossed (shaded areas; data for A/Ec = 5). This leads to the 
formation of metallic puddles, i.e., extended regions which are 
either n- or p-doped. An example is shown on the right (T = 0, 
A/Ac = 10, green/brown: n/p-doped). With increasing tem¬ 
perature, the fluctuations of the potential decrease (dashed 
lines in left panel) due to screening by thermally activated 
carriers thereby suppressing puddle formation. 


BiSbTeSe 2 and the sensitivity of transmittance measure¬ 
ments to weak absorption features. The corresponding 
spectral weight is strongly temperature dependent at low 
temperatures. Based on numerical simulations, we will 
argue that this temperature dependence is indeed char¬ 
acteristic of the mechanism of puddle formation by fluc¬ 
tuations of the Coulomb potential. 


II. EXPERIMENTAL RESULTS 
A. Optical spectroscopy 


The complex optical conductivity ai{uj) + iu 2 {io) of 
single-crystalline BiSbTeSe 2 was determined from in¬ 
frared transmittance and reflectance data which were 
complemented by ellipsometric measurements at higher 
energies (see Methods, Sec. VIB). An overview of cti (w) 
in the infrared range is plotted in Fig. on a logarith¬ 
mic scale. The spectra reveal the steep increase of cri(a;) 
caused by the onset of excitations across the gap A. At 
5K, we find A = 0.26eV (2100cm"!). At 300K, A is re¬ 
duced by about 40%, it decreases with a slope of roughly 
3.6cm“^/K. Similar results for the temperature depen¬ 
dence were reported for related topological insulators. 
For more details, see Supplemental Material [53] . 


The main focus of the present study is, however, on the 
electronic contribution to the optical conductivity below 
the gap and its peculiar temperature dependence. In 
the temperature window from 40 to 60 K, CTi(a;) reaches 
values as low as 0.3 (Hcm)“^. Most remarkably, the tem¬ 
perature dependence of ai{u}) is highly non-monotonic. 
In the frequency range of about 300-1100cm“^, cri(w) 
is more than three times larger at 5K than at 50 K. The 
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FIG. 2: Optical conductivity of BiSbTeSe2 on a logarith¬ 
mic scale. Weak absorption features below the gap with 
0-1(0;) < 10(flcm)“^ were obtained from the transmittance 
for a sample thickness of d = 102 ytm, while data with 01 (oj) > 
20 (flcm)“^ in the opaque range were derived via a Kramers- 
Kronig analysis of the reflectivity. In combination, these data 
sets give an excellent account of oi(a;). 


d = 102, 130, and 183 /xm agree very well with each other 
within the experimental uncertainty (see Fig. S3 in Sup¬ 
plemental Material [53]). This proves the bulk character 
of the excitations in the investigated frequency range. 
Theoretically, one may expect two contributions from the 
surface state: a Drude peak arising from surface conduc¬ 
tion and interband excitations within the Dirac bands. In 
BiSbTeSe2, the Fermi level is close to the Dirac point |25j . 
giving rise to a small density of surface states. Moreover, 
Dirac fermions show a large mobility. The respective 
narrow Drude peak is located below the frequency range 
addressed in our data, in agreement with terahertz data 
on thin films of 612803 and Bii.aSbo.sTei,3861.2 [SHI El]. 
Interband excitations within the Dirac bands contribute 
at higher frequencies. For huj> 2Ep, a universal conduc¬ 
tance Go = 7re^/8/i « 1.5 • 10“®/D has been predicted 
[m [14]. For d=100^m, this is equivalent to a bulk 
conductivity of 0.0015 (Dcm)“^, which is two orders of 
magnitude smaller than the lowest values observed in 
Bi8bTe8e2, see Fig. We therefore conclude that all 
of our observations reflect bulk properties. 


C. Electronic contribution to a\(uj) 


rise of (Ti (w) upon heating above 50 K agrees with the DC 
conductivity a 1(10 = 0) measured in transport [13], but 
the increase of cti(w) upon cooling below 50 K strongly 
deviates from the transport results. This discrepancy of 
DC and optical conductivities is the smoking gun for the 
puddles, as it is a natural consequence of the carrier lo¬ 
calization within puddles. 

Note that for all temperatures the measured values of 
tTi(a;) below the gap are by far the lowest reported thus 
far for the entire family of Bi2_x8bxTe3_ySej,. In Bi2Te3 
and 612803, the Drude contribution of extrinsic carri¬ 
ers with typical densities iV«10^®cm“^ extrapolates to 
DC values of cri(O) « 1000 (Dcm)“^ [THl |H1133 EHl 
El]. In compounds with smaller N such as Bi2Te2Se, 
impurity absorption bands with peak values of 50- 
100(Dcm)“^ were reported [^SHSO] , one to two orders 
of magnitude larger than the conductivity observed by 
us. Such pronounced impurity bands are apparently ab¬ 
sent in BiSbTeSe2, in agreement with recent reflectiv¬ 
ity data |52] . which were, however, not sensitive enough 
to reveal the comparably weak absorption features with 
ai{uj) < 10 (Dcm)“^ observed by us in transmittance. 


B. Absence of surface contributions 


An important question is whether the spectral weight 
observed below the gap can be related to the surface 
states of the topological insulator. This can, however, 
be excluded by comparing data for different thicknesses 
d obtained successively on the same sample (see Meth¬ 
ods, Sec. VI). At each temperature, results for cri(a;) for 


Figure shows cri(w) on a linear scale for frequencies 
below the gap. Several contributions can be identified in 
this frequency range, see inset of Fig. [^ Below about 
150cm“^, (Ji{ui) is dominated by a phonon contribution 
with a peak value of the order of 10^ (Dcm)“^ [501 [52] 
which can be identified in the reflectivity data (see Fig. 84 
in Supplemental Material [33]). Above 150cm“^, we find 
a tiny absorption band extending up to about 350 cm“^ 
with a peak value of about l(Dcm)“^. Based on the 
frequency range and the tiny spectral weight, this can 
be attributed to a multi-phonon contribution, i.e., two- 
and three-phonon excitations. The remaining contribu¬ 
tions of electronic origin we fit with a tiny, temperature- 
independent constant term of about 0.2(Dcm)“^ and a 
strongly temperature-dependent Drude peak. Well above 
50 K, the interpretation of this feature as a Drude peak of 
thermally activated carriers is supported by the absolute 
value of cri{ijj), by the peak width, and by the tempera¬ 
ture dependence of the spectral weight, as shown below. 
The main focus of our study is, however, on the reap¬ 
pearance of spectral weight at low temperatures, which 
can be attributed to locally n- or p-doped puddles. The 
optical conductivity of such puddles is also expected to 
be of Drude form for frequencies above a cut-off lUc given 
by the Thouless energy, determined by the time scale 
needed to diffuse through a puddle. Due to the large 
size of the puddles, the cut-off uic is orders of magnitude 
smaller than the frequency range investigated by us. Ac¬ 
cordingly, we fit the data using the Drude model also at 
low temperatures. 

In the Drude model, cri(a;) depends on the scat¬ 
tering rate 1/r and the effective carrier density 
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compound 

Aeff[10i Vernal 

l/r[10^Vs] 

T[K] 

Ref. 

BiSbTeSe 2 

0.02/0.6 

140 

50 / 300 

this work 

Bi2Te2Se 

1.9 

40 

300 

[48] 

Bi2Se3 

2.9; 18 

4; 23 

6; 300 

m-: [55] 

Bi2Te3 

33; 46 

4.7; 5.6 

10; 10 

m -: [SSI 


TABLE I: Effective carrier densities NeS = Nnie/m* and 
scattering rates for different compounds. Carrier densities 
from Refs. | 18[ 148 1 were calculated from the unscreened 
plasma frequencies given there. In Ref. [5S], the screened 
plasma frequency ujp/y/W^ is given together with Soo = 29.5 
for Bi 2 Se 3 . 


Neff = Nme/m*, 


ai (w) 


c^i(O) 
1 + 


A^eff rime 
1+Uj‘^T^ ’ 


( 1 ) 


where e and rUe denote charge and mass of a free elec¬ 
tron, respectively, and m* is the effective band mass. 
Well above 50 K, the fit results for cri(O) are consistent 
with DC resistivity data of samples with the same stoi¬ 
chiometry |24) . 

Comparing our result for Neff at room temperature 
with Hall-effect data [24], we find m* jme =0.2, see Sup¬ 
plemental Material |53j . This agrees with results for 
Bi 2 Se 3 , where values between 0.14 to 0.24 were de¬ 
rived from the cyclotron mass of the bulk conduction 
band depending on the orientation of the cyclotron orbit 
PP] . Using m*/me = 0.2, we deduce a carrier density as 
low as iV « 4 • 10^® cm“^ between 40 K and 60 K. The 
temperature-driven increase of Neff above about 50 K 
can be described as activated behavior with an activa¬ 
tion energy i?^ = 26meV « 0.1 A, see inset of Fig. 
This agrees with Ea = 22-30 meV derived from trans¬ 
port measurements on BiSbTeSe 2 for temperatures above 
100 K [24]. The small activation energy has been pro¬ 
posed to be a clear signature of strong Coulomb fluctua¬ 
tions [5P1157] . 

From the peak width we obtain 1/r « 1.4 • 10^^ s“^ 
roughly independent of temperature as expected for a 
scattering mechanism arising from the random posi¬ 
tion of defects. With m* jme = 0.2, this corresponds 
to a mobility p = eT/m* « 70cm^/Vs, in excellent 
agreement with the value of 73cm^/Vs from Hall data 
[21] on Bii.sSbo.sTei, 7801 . 3 . Scattering rates of differ¬ 
ent compounds are compared in Table |l] Compensated 
BiSbTeSe 2 shows the smallest carrier density and by far 
the largest value of 1 /r which supports that defect scat¬ 
tering is dominant. 


D. Puddles 

Our main result is the dramatic reappearance of low- 
frequency spectral weight at temperatures below 50 K, 
see Fig. [^ The charge carriers responsible for this do, 
however, not contribute to the DC conductivity, and 



Frequency [cm'] 


FIG. 3: Optical conductivity below the gap. At 40-50K, 
o\{ijj) is tiny below the gap. With increasing tempera¬ 
ture, we identify a Drude peak of activated carriers with a 
strongly temperature-dependent spectral weight (see bottom 
panel and Fig.|^ and a large and approximately temperature- 
independent scattering rate, 1/r ~ 1.4 • lO'^'^s”^. Most re¬ 
markable is the reappearance of low-frequency spectral weight 
below about 50 K, which reveals the formation of puddles, see 
top panel. The inset shows a fit of the 40 K data with four 
contributions: a phonon at 70cm“^, a multi-phonon band 
at 275cm“^, a constant background of 0.23(Dcm)“^, and a 
broad low-frequency band for carriers localized within puddles 
(red). The phonon and the background are kept constant in 
the fits of other temperatures (dashed lines in both panels). 


cri(a;) at 5K is about an order of magnitude larger than 
cji{uj = 0) [24]. This is consistent with a picture of well 
separated metallic puddles contributing to cri(w) but not 
directly to DC transport. 

The effective carrier density amounts to A/ff,p « 1.2 • 
10^®cm“^ at 5 K. Using the value of the effective mass de¬ 
termined at 300 K, this corresponds to an average carrier 
density N-p « 2 • 10 ^^cm“^, which is, however, expected 
to be distributed in a highly non-uniform way due to 
puddle formation. With increasing temperature, the car- 
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FIG. 4: Effective carrier density Symbols depict fitting 
results for the low-frequency absorption band obtained for 
different sample thicknesses d. Below about 50 K, the carriers 
can be attributed to puddles. Solid line: activated behavior 
with an activation energy Ea ~ 26meV. Inset: same data on 
a log scale vs. 1/T. 


rier density shows a rapid drop by a factor of 4-6 at 
a temperature scale of the order of 30-40K, see Fig. 
Below we will show that this temperature scale has to be 
identified with the energy scale which agrees quanti¬ 
tatively with theoretical expectations. Also the average 
carrier density Np can be explained by our numerical 
simulations. 

Note that an unconventional - but much weaker - tem¬ 
perature dependence of the carrier density has been ob¬ 
served before in this family of topological insulators. In 
Bi 2 Te 3 , an unconventional decrease of A^efi by up to a 
factor of 2 has been observed between 5K and 300 K 
[Ml 1551 157] . For compensated Bi 2 Te 2 Se with Neff ~ 
10 ^® cm“^, a non-monotonic behavior of A^eff with a min¬ 
imum in the range of 50 K to 150 K was reported ISD- 
However, Neff in Bi 2 Te 2 Se changes by less than 10% be¬ 
tween 5 K and 50 K and by about 20 % between 50 K and 
300 K, whereas we find a drastic change by more than a 
factor of 10 in BiSbTeSe 2 , see Fig. We emphasize that 
our results are based on samples with very low carrier 
density in combination with the enhanced sensitivity for 
weak absorption features offered by transmittance mea¬ 
surements. 


III. MODELLING THE FORMATION AND 
DESTRUCTION OF PUDDLES 


Following Skinner, Chen, and Shklovskii [55], we use 
a simple classical electrostatic model to describe the 
formation of puddles in a compensated semiconductor. 
The model assumes that donors and acceptors are lo¬ 
cated at random positions in space. Their average 
densities are given by No and Na, respectively, with 


Ndef = {Na + Njy)/2. We are mainly interested in the 
experimentally relevant limit of almost perfect compen¬ 
sation where K = Na/Njj is close to 1. The binding en¬ 
ergy of charges to defects is small due to the large dielec¬ 
tric constant (see Fig. S4 in Supplemental Material [55]), 
thus donors and acceptors are shallow with energy levels 
very close to ±A/2. This situation is described by the 
Hamiltonian 

H — — fiUi + - Vr--rjqiqj ( 2 ) 

i hj 


where Ui = 0,1 denotes the number of electrons on a 
donor (/i = l) or acceptor (/j=-l) site. The charge of 
a donor (acceptor) amounts to qi = l {qi=-l) if it has 
donated (accepted) an electron to (from) another de¬ 
fect, otherwise defects are charge neutral, qi = 0. The 
Coulomb potential is supplemented by a short-distance 
cutoff qb, Vri-rj = e^/{ 47 r£o£(|ri — -|- a^)^/^}, which 

effectively takes into account the finite extent of the 
wave functions of the shallow impurity states [36] . The 
value of Os turns out to have little influence [55] and is 
set to a_B = (2/Adef)^^^ for all of our simulations. Ex¬ 
pressing all distances in units of the average distance 
of dopants, ddef = l/-^def > energies in units of 

the Coulomb interaction between neighbouring dopants, 
Ec = e^/(47r£o£ddef), all properties of the model depend 
on A/Ec, K, and T/Ec- 

The strength of Coulomb interactions and therefore 
Ec strongly depend on the dielectric constant e which is 
strongly frequency dependent in BiSbTeSe 2 and related 
compounds. Below the gap but above the phonons, we 
find £«35, which increases to £ « 200 for a; —>■ 0 due to 
a huge phonon contribution, see Fig. S4 in Supplemental 
Material [55]. As puddles are static objects around which 
the highly polarizable ions will adjust their positions, the 
a; —)■ 0 value £ « 200 should be most relevant for our 
model and is therefore used in the following. 

Besides the donor and acceptor states, no further con¬ 
duction or valence electron states are taken into account 
in Eq. ([^. For a gap of A/fcs^OOOOK, the contribution 
of intrinsic carriers thermally activated across the gap can 
be neglected at low temperatures. Also the intrinsic car¬ 
rier density within a puddle can be neglected. This is due 
to the small effective mass m* in combination with the 
small value of the Fermi energy Ep^Ec within the pud¬ 
dles (see below). Using m*lme = 0.2 and Ep = 50K for 
a single spherical band, one obtains an electron density 
of 10 ^^cm“^, more than an order of magnitude smaller 
than the typical density of defects. 

While the classical model of Eq. ([^ is strongly sim¬ 
plified, it is a powerful tool [55H55] to obtain a semi- 
quantitative understanding of puddle formation at T = 0. 
We will show below, that it also describes the destruction 
of puddles with increasing temperature. Most impor¬ 
tantly, the model is sufficiently simple to allow for quan¬ 
titative numerical simulations both at T = 0 and at finite 
T (see Methods, Sec. VIC). We are able to obtain results 
with only small finite-size effects for values of A/E^ up 
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FIG. 6: Puddle formation and doping. At T = 0, the fraction 
Pp = E„> 4 Po(n) (see Fig. Ef of dopants organized in puddles 
rises strongly in systems with small deviations |1 — A”! ^ 1 
from perfect compensation. The scaling plot of ppA/Ec as 
function of (1 — K) A/Ec = for different values of 

A/Ea shows that pp « 0.06 Ac/A for perfect compensation, 
while Pp ~ 0.33 (1 — K) for 1 — A > 0.3 Ac/A. 


FIG. 5: Identification of puddles. Top panel: Schematic pic¬ 
ture of an electron puddle in a region where the density of 
donors (green circles) is larger than the density of acceptors 
(yellow circles). The symbols -f, —, and 0 indicate the dopant 
charge. To identify a puddle, we consider a sphere of radius 
ro = 1.42 ddef around each neutral dopant. Then, we count 
the number n of neutral dopants of the same kind within 
each sphere to obtain the distribution function po{n) (normal¬ 
ized by the total number of dopants). This is shown in the 
lower panel for A/Ac = 15 (left: perfect compensation, A = 1; 
right: A = 0.95) and various temperatures. We identify the 
fraction Pp of dopants located well within a puddle with 
Pp — En> 4 Po(n), i.e., considering all neutral dopants with 
at least four neutral neighbours (cf. vertical dashed lines). At 
r = 0 (solid blue line), most neutral dopants are organized in 
puddles, having a substantial number of neutral neighbours. 
With increasing T, the number of neutral dopants with many 
neutral neighbours drops for T < Ec, which is a clear sign 
for the destruction of puddles (see Fig. [^. The number of 
isolated neutral dopants with no neutral neighbours rises in¬ 
stead. 



FIG. 7: Destruction of puddles with increasing temperature. 
The fraction pp (see Fig. of dopants organized in puddles 
drops rapidly as function of temperature at a temperature 
scale set by the Goulomb interaction Ac between neighbour¬ 
ing dopants. Numerical results are given for A/Ac = 9, 12, 
and 15 (left: perfect compensation, A = 1; right: A = 0.95). 
The scaling plots {ppA/Ec for A = 1 and Pp/{1 — A) for 
1 — A > 0.3 Ac/A, see Fig. as function of TjE^ allow us 
to extrapolate the results to larger values of A/Ac. Scaling 
demonstrates that the destruction of puddles always occurs 
at approximately the same value of A/Ac. 


to 15, see Fig. S5 in Supplemental Material [53]. Scaling 
arguments then allow us to address the experimentally 
relevant regime of A/Ac < 100 (see Sec. IV). 

At r = 0, we reproduce the results of Ref. [35], mini¬ 
mizing the energy by a pairwise exchange of charges. We 
extend, however, the simulation to finite temperatures 
using a Monte Carlo approach (see, e.g.. Ref. [60] for fi¬ 
nite T simulations for other Coulomb systems). Puddles 
are formed from occupied donor states or empty acceptor 
states, see Fig. [^ These correspond to neutral dopants, 
where, e.g., an electron compensates a positively charged 
donor ion. To detect puddles in our simulation, we distin¬ 
guish extended regions of neutral dopants from isolated 
sites. Around each neutral dopant, we consider a sphere 
of radius rg and count the number n of other neutral 
dopants of the same type, see the sketch in Fig. [^ The 
radius rg = 1.42 duef is chosen such that on average there 


are 12 dopants (the number of nearest neighbours for 
close-packed spheres) of the same type within the sphere. 
With po{n) we denote the fraction of dopants which are 
neutral and have n neutral neighbours. This fraction 
is plotted in Fig. for perfect compensation A = 1 (left 
panel) as well as for A = 0.95 (right). For r = 0, most 
neutral dopants are organized in puddles, i.e., have many 
neutral neighbours. Note, however, that Poin) at T = 0 
is peaked at a value of n substantially smaller than 12 
which implies that an impurity band formed by donor 
or acceptor states within a puddle is only partially filled. 
This is related to the fact that the energy scale governing 
puddle formation (the depth of the potential A+(r) — p, 
in Fig. is given by and is of the same size as the 
fluctuations of energy of neighbouring dopants. 

Upon increasing T, the total density of neutral dopants 
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{Na + Nd) J2‘^=oPo(.^) increases. For low T this in¬ 
crease is due to the thermal activation of carriers outside 
of the puddles, as Po{n) rises only for small number n 
of neighbours, see Fig. At the same time, the num¬ 
ber of neutral dopants with large n decreases, thus also 
the number of carriers inside the puddles decreases. The 
reason is that the thermally activated charges screen the 
Coulomb potential, which leads to a pronounced reduc¬ 
tion of the fluctuations of the Coulomb potential, see 
Fig-lU and therefore to the destruction of puddles in a 
highly non-linear process. As we will show, this mecha¬ 
nism of puddle destruction is remarkably robust against 
changes of IS./E^. and deviations from perfect compensa¬ 
tion, K = 1. The identihcation of this mechanism and of 
the corresponding energy scales is the main result of our 
theoretical analysis. 

To quantify this effect, we count the fraction pp = 
Sn>no 7 'o(’^) of neutral dopants located well within a 
puddle. For the discussion, we choose no = 4, i.e., we 
count those neutral dopants with at least four neutral 
neighbours, as indicated by the vertical dashed lines in 
Fig. We have also used other values of uq (not shown) 
to confirm that the results do not depend qualitatively 
on this choice. A similar approach to identify clusters in 
random systems is, for example, used in Ref. [6T] . 

We first consider the limit T —>■ 0. A scaling collapse 
of all results is obtained when pp ■ A/Ec is plotted as 
function of (1 — K)A/Ec, see Fig. For K = l, the 
density of neutral dopants contributing to puddles is of 
the order N^gfEc/A with a small prefactor. This den¬ 
sity does, however, rapidly increase when small devia¬ 
tions from perfect compensation are taken into account. 
This occurs roughly when \N£) — Na\, the uncompensated 
part of the doping, becomes larger than the density of 
neutral dopants organized in puddles at K = 1. The lin¬ 
ear growth of Pp with |1 — iF| corresponds to the effect 
that adding, e.g., a small amount of extra donors to the 
perfectly compensated system does not give rise to a uni¬ 
form doping but instead increases the number of neutral 
donors organized in puddles. This non-uniform doping 
originates again in the large-scale inhomogeneities of the 
Coulomb potential and is important also for the behavior 
at finite temperature. 

Upon increasing T, the fraction pp of dopants in pud¬ 
dles drops sharply, as shown in Fig. For perfect com¬ 
pensation, plots of Pp{T) ■ A/Ec versus ksT/Ec are in¬ 
dependent of A/Ec for T,Ec A (left panel of Fig. [^. 
This allows us to predict that the Coulomb interaction 
Ec between neighbouring dopants sets the temperature 
scale for the destruction of puddles even in the experi¬ 
mentally relevant regime A/E^ ~ 100 (see below). The 
right panel of Fig. where Pp/{1 — K) is investigated 
for K = 0.95, shows that this physics is not affected by 
small deviations from perfect compensation. In the in¬ 
vestigated parameter regime (|1 — AT] ^ 1 and A Ac), 
we find that the temperature scale for the destruction of 
puddles is independent of doping effects and always given 
by Ec- This has to be contrasted with the strong effects 


of doping on the absolute value of Pp which changes by 
an order of magnitude, see Fig. 

While the destruction of puddles with increasing tem¬ 
perature has, to our knowledge, not been investigated 
before, our results are fully consistent with known prop¬ 
erties [36l |39] of such Coulomb systems. For T = 0 and 
AT = 1 , it was observed [SHIES] that the density of states 
of effective single-particle levels is roughly constant (up 
to the famous Efros-Shklovskii Coulomb gap at low fre¬ 
quencies) in an energy window set by ± (^ + F'c) and 
therefore given by 1 /A for A ^ E^- Furthermore, the 
energy scale of the carriers within a puddle is set by 
Ec- In combination, this implies that the fraction of 
charge carriers in puddles is of the order of E^j A in the 
case of perfect compensation, as observed in our sim¬ 
ulations. This value is, nevertheless, surprisingly large 
when compared to the much smaller fraction of charges, 
~ (Ac/A)^, which should be sufficient to compensate 

charge fluctuations of order ^A/jef Ag within the non¬ 
linear screening radius Rg discussed in the introduction. 

IV. DISCUSSION 

Both of our main experimental observations, the pres¬ 
ence of a sizable optical weight at low T and its rapid 
drop on a small temperature scale of the order of 30- 
40 K, agree qualitatively with our numerical simulations 
based on a model of shallow donors and acceptors in¬ 
teracting by long-ranged Coulomb interactions. The re¬ 
maining task is to compare the parameters of theory and 
experiment quantitatively. 

For the Coulomb energy between neighbouring 
dopants, theory predicts 
2 

Ec = • 20 - 40 K , (3) 

where we used e « 200 (see Fig. S4 in Supplemental Mate¬ 
rial [53] ) and assumed that the density of shallow donors 
and acceptors is in the range of Vdef ~ 10 ^® — 10 ^° cm“^ 
as estimated from the carrier density in uncompensated 
samples, see Tab. I. This is fully consistent with the 
experimentally observed temperature scale of 30-40K, 
which translates to a ratio of A/Ec « 75 — 100. This 
agreement does not only corroborate the assertion that 
puddles exist in BiSbTeSe 2 but also confirms the scenario 
that they arise from strong fluctuations of the Coulomb 
potential. Moreover, this agreement indicates that the 
optical determination of Ec may turn out to be a use¬ 
ful tool to estimate the defect density Vdef, a quantity 
which is difficult to assess in a compensated semiconduc¬ 
tor. This interesting result will have to be tested in future 
experiments. 

A quantitative prediction for the effective carrier den¬ 
sity NcS is more subtle as this quantity depends for 
large A/Ec sensitively on the precise amount of com¬ 
pensation which is not known experimentally. For per¬ 
fect compensation, Ar = l, our results for pp show that 




Np/Ndei should take values of the order of 0.1 Ed A. ~ 
10“^. The defect density A^def can be estimated from 
the experimental value for Ec and Eq. @ which gives 
.^def = 5 —10 • 10^® cm“^. Combining this with the exper¬ 
imental estimate Np « 2.4-10^^ cm“^ at 5K (see above), 
we obtain Np/N^ei ~ 0.002 — 0.005, a factor 2—10 larger 
than the theoretical estimate for itT = 1. Most likely, this 
just means that compensation is not perfect in our sam¬ 
ple. As can be seen from Fig. tiny deviations from 
perfect compensation at a level of 1 — AT < 1.5Ec/A 
can easily explain the observed spectral weight. With 
AjEc « 75 — 100, a deviation from perfect compensa¬ 
tion of just 1 or 2 % is sufficient to obtain a consistent 
description of the experiment. 

Taking both the extremely simplified nature of the 
theoretical description and the uncertainties in parame¬ 
ters like the effective mass into account, the quantitative 
determination of the parameters should perhaps not be 
taken too literally. They are, however, highly plausible, 
suggesting that at low temperatures we achieve agree¬ 
ment between theory and experiment at least on a semi- 
quantitative level. 


V. CONCLUSION 

Our optical conductivity data of the almost perfectly 
compensated topological insulator BiSbTeSe 2 reveal the 
existence of puddles at low temperatures as well as their 
destruction on a temperature scale of 30-40K. Both the 
spectral weight and the temperature scale agree semi- 
quantitatively with our numerical simulations based on 
a model of shallow donors and acceptors interacting by 
long-ranged Coulomb interactions. We have shown that 
puddles are suppressed by thermally activated charges 
which screen the Coulomb potential. The temperature 
scale of puddle destruction is set by the Coulomb interac¬ 
tion Ec between neighbouring dopants. This mechanism 
works both for near-perfect and perfect compensation. 

Puddle formation driven by long-ranged Coulomb in¬ 
teractions is not only of importance in compensated semi¬ 
conductors but also for other materials with a vanish¬ 
ing density of electronic states including Dirac matter 
in two or three dimensions, like graphene |62] or Weyl 
semimetals. For the physics of topological insulators, 
puddle formation in compensated samples has both pos¬ 
itive and negative effects. While the strong fluctuations 
of the Coulomb potential imply that it is more difficult 
to reach high bulk resistivities despite of perfect compen¬ 
sation, they can also help to localize electrons or holes in 
puddles in situations where the compensation of donors 
and acceptors is not perfect. The surface states of topo¬ 
logical insulators can provide extra screening channels, 
thus suppressing puddle formation close to the surface 
or in thin samples. It will therefore be interesting to 
study both experimentally and theoretically, how puddle 
formation depends on sample thickness and other param¬ 
eters and how it interacts with the charge density of the 


topological surface states. Controlling puddle formation 
may turn out to be a key step for further reduction of 
bulk transport In topological insulators. 


VI. METHODS 
A. Samples 

The compound BiSbTeSe 2 belongs to the family 
of A 2 B 3 tetradymites (A = Bi,Sb; i? = Te,Se) showing 
rhombohedral structure (space group R3m) [Ml (53] with 
three quintuple layers per unit cell stacked along the [ 111 ] 
direction. Single crystals of BiSbTeSe 2 were grown start¬ 
ing from high-purity elements as described in Ref. [24j . 
The crystals were cut into platelets with typical dimen¬ 
sions of 3 X 3mm^ within the (111) plane. Due to the 
weak van der Waals bonding between quintuple layers, 
the samples can be cleaved easily along the ( 111 ) plane 
using adhesive tape. This yields shiny plane-parallel sur¬ 
faces. 

The defect density can be reduced by chalcogen or¬ 
der - A- - A- B^ within the quintuple layers m 
|23||24l Eg. In Bl 2 Te 2 Se with Te-Bi-Se-Bi-Te order 
[niEg, the B^^ sites are exclusively occupied by Se ions. 
In the solid solutions Bi 2 _a;Sba;Te 3 _ySej, with y > 1, 
the composition has been optimized with the aim to 
achieve full compensation |24l EH]- Chalcogen order is 
preserved to some extent, as shown by x-ray diffraction 
[Mj . Among these solid solutions, BiSbTeSe 2 was re¬ 
ported to show the highest DC resistivity at 2 K [M] . It 
reaches 3 Dcm but varies by about a factor of 3 for sam¬ 
ples with the same nominal composition. This can par¬ 
tially be attributed to a thickness dependence [M| EZl [HH] 
related to a finite conductance contribution of the surface 
but also reflects different defect concentrations [M] . 


B. Optical measurements 

Infrared reflectance and transmittance measurements 
were performed with unpolarized light in the frequency 
range of 50-7500cm“^ (6meV-0.93eV) using a Bruker 
IFS 66 v/S Fourier-transform spectrometer equipped with 
a continuous-flow He cryostat. The transmittance 
was recorded at normal incidence with the electric field 
parallel to the cleavage plane, while the reflectivity RijE) 
was measured under near-normal incidence. Addition¬ 
ally, ellipsometric data were obtained using a rotating 
analyzer ellipsometer (Woollam VASE) equipped with a 
retarder between polarizer and sample. The ellipsomet¬ 
ric data were collected at room temperature in the pho¬ 
ton energy range of 0.75-5.5eV (6050-44360cm“^) for 
three different angles of incidence (60°, 70°, and 80°). 
Reflectance data and ellipsometric data were measured 
on a sample with a thickness of dw 1.1 mm. For the trans¬ 
mittance measurements, we started on a sample with 
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a thickness of d = (183 ± 5)/^m (see Fig. SI in Supple¬ 
mental Material [53]). The value of d was determined 
mechanically using a micrometer screw. For this rather 
thick sample, the accuracy of 5 /rm corresponds to an er¬ 
ror of 2.7%. Subsequently this sample was cleaved sev¬ 
eral times using adhesive tape, and the transmittance 
was measured successively on the same sample for a se¬ 
ries of different thicknesses, d = 183, 130, and 102 /im. 
The latter two values were determined by comparing the 
Fabry-Perot interference fringes which arise in a trans¬ 
parent frequency range due to multiple reflections within 
the sample. Due to the shiny and plane-parallel surfaces 
obtained by cleaving, the interference fringes are partic¬ 
ularly pronounced, see Fig. SI in Supplemental Material 
[55] . For two samples a and h with different thicknesses 
da and df,, the thickness ratio can be determined from 
the fringe periods, ^Vaj = db/da, with an accuracy of 
better than 0.5%. This is important for the comparison 
of results obtained for different thicknesses, see Fig. S3 
in Supplemental Material [53], and thus for the question 
whether there is a finite contribution of surface states. 

In the transparent frequency range, the complex opti¬ 
cal conductivity d(uj) = cfi{uj) -I-* 0 - 2 ( 0 ;) was determined 
from T(uj) and R{ui) [^. In the opaque range, dluj) was 
obtained via a Kramers-Kronig analysis of R{u}), which 
at high frequencies was extrapolated using the ellipso- 
metric results. 


C. Monte Carlo simulations 

The model defined in Eq. ([^ is simulated at finite T 
with a standard Monte Carlo algorithm (Metropolis). 
Periodic boundary conditions for the Coulomb poten¬ 
tial are imposed by using always the shortest distance 
on the 3—torus for its computation. We start the sim¬ 
ulations from a configuration where all Na acceptors 
and (1 — K)Njj donors are occupied, such that the to¬ 
tal system is charge neutral. We only consider config¬ 


urations which keep charge neutrality by using a pair¬ 
wise exchange of charge in each Metropolis step. For 
T —>■ 0 we average over 100 disorder realizations. At 
high temperatures, averages over only 10 different real¬ 
izations turn out to be sufficient. We simulate up to 
Na + Nu = 2iVdef = 2 X 38^ « 110.000 dopants. For the 
parameters K = l, T/Ea = 0, and A/Ea = 15 finite-size 
effects are largest, see Supplemental Material [53] for de¬ 
tails. Therefore, we used 2 x 60^ = 432.000 dopants 
for this particular set of parameters (a single triangle in 
Figs. and [^. 

For T = 0 the algorithm is identical to the one used 
in Ref. [35]. It yields only a local and not a global 
minimum of the (free) energy but such pseudo ground 
states are known to describe the properties of real ground 
states with high accuracy |36j . In contrast to simula¬ 
tions with local interactions, the numerical costs for the 
long-ranged Coulomb interactions increase strongly for 
increasing temperature and scale with {NaN: each 
update of the charge configuration implies that the local 
energies of all other dopants have to be recomputed (see 
Ref. [36]). As puddle formation occurs at length scales of 
(A/£’c)^/A^def, the number of dopants needed in a simu¬ 
lation grows with (A/Ac)® and a worst-case estimate for 
the computational costs is (A/i?c)^^ for finite tempera¬ 
tures. 
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Supplemental Material: 

Revealing puddles of electrons and holes in com¬ 
pensated topological insulators 


band-structure calculations for quaternary compounds 
have not been reported thus far. 


I. TRANSMITTANCE SPECTRA 


III. OPTICAL CONDUCTIVITY FOR 
DIFFERENT SAMPLE THICKNESSES 


Figure SI shows the transmittance T{ui) measured on 
a sample with a thickness of d = 183 /rm. At 300 K, 
exceeds the noise level only between 830 and 1335 cm“^ 
and stays below 0.5%. However, T{w) strongly increases 
upon cooling down. At 50 K, the sample is transparent 
between about 180 cm“^ and 2060 cm“^. Most remark¬ 
ably, T[uj) decreases upon further cooling below about 
50 K. At 5K, T{uj) is much lower than at 50 K for fre¬ 
quencies not too close to the gap. In the highly transpar¬ 
ent range, T(u}) exhibits Fabry-Perot interference fringes 
(see right panel of Fig. 


between front and bac 


SI I caused by multiple reflections 


:k surface. The fringes are partic¬ 
ularly pronounced due the shiny, plane-parallel surfaces 
of cleaved samples. 


II. TEMPERATURE DEPENDENCE OF THE 
GAP 


Figure S2 depicts A(T) as determined from the onset 
of transmittance for a sample thickness of d = 102 /im. 
The gap shifts by almost 40% between 5 K and 
300K, from about 2112cm“^ (262meV) to 1292 cm“^ 
(160meV). A ht based on the empirical Varshni equa¬ 
tion A(T) = A(0) — /3 • r^/(T -I- To) (red line) yields 
A(0)=2125cm-\ To = 86K, and ^ = 3.6cm-VK or 
5.2 ke- Our result for (3 is framed by values reported 
for binary compounds, e.g., 0.77cm“^/K for Bi 2 Te 3 [HD, 
1.6cm ^/K and 2.0cm ^/K [S5] for Bi 2 Se 3 , and 
5.6cm“^/K for Sb 2 Se 3 [S2]. Band-structure calculations 
for binary compounds indicate that A is very sensitive 
to temperature [S4l[S5j. To the best of our knowledge. 


The bulk conductivity (Ji{u}) and a possible surface 
conductance G(l>j) can be disentangled by comparing re¬ 
sults for different thicknesses d. If the data analysis is 
based on the assumption of a bulk-only character, a fi¬ 
nite surface conductance gives rise to an apparent thick¬ 
ness dependence of (Ti(w), i.e., the calculated cri(a;) in¬ 
creases upon decreasing d. In contrast, our results for 
ai(u}) for different d agree very well with each other, see 


Fig. 

0.99 


S3 At 5 K, the ratio 


lies within the range 


1.05 below 1500 cm“^. This strongly supports a 
bulk-only character of the investigated excitations. At 
1000cm“^, the difference between data for d=183/xm 
and 102 /rm lies in the range —0.1 Bcm“^ to -1-0.2 f2cm“^ 
for different temperatures. This can be considered as 
the experimental uncertainty, which can be attributed to 
small experimental errors concerning the absolute value 
of T(tLi) or the thickness d. In our data, fTi(a;) typically 
takes the smallest value for d = 130 /rm. These values are 
in particular smaller than the results for d = 183 /rm, in 
contrast to the expectations. This points towards a small 
systematic error of the absolute value. 


IV. EFFECTIVE MASS 

The effective mass m* can be estimated by comparing 
our results with Hall-effect data |S6) which yield a Hall 
constant of Rh = 5cm^/C at 300 K. In a perfectly com¬ 
pensated semiconductor, both electrons and holes con¬ 
tribute to transport. In this case, Rh gives an upper esti¬ 
mate of the carrier density, V = 1.25 • 10^® cm“® at 300 K. 
However, our analysis (see Discussion in main text) sug- 
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FIG. SI: Infrared transmittance spectra of BiSbTeSe 2 . Left: 
T{u}) measured on a sample with a thickness of d=183pm. 
Right: same data on an enlarged scale, highlighting the pro¬ 
nounced interference fringes. 


FIG. S2: Temperature dependence of the gap. Symbols depict 
A(r) as determined from the onset of transmittance for a 
sample thickness of d = 102 fim. The red line shows a fit based 
on the empirical Varshni equation. 
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Frequency [cm '] 


FIG. S3: Results for the optical conductivity cri{u}) for dif¬ 
ferent sample thicknesses d agree very well with each other, 
strongly supporting a bulk-only character. 


gests that the sample of BiSbTeSe 2 shows small devia¬ 
tions of 1-2% from perfect compensation, i.e., one car¬ 
rier type predominates. We thus expect that the result 
for N derived from Rh is reliable. From the optical data, 
we find an effective carrier density ■ 10^® cm“® at 

300K, from which we obtain m* In BiSbTeSe 2 , 

the cyclotron mass of the bulk bands is unknown since the 
measured quantum oscillations arise only from the sur¬ 
face states [SB]. However, the cyclotron mass of the bulk 
conduction band has been measured in Bi 2 Se 3 , showing 
m*/me between 0.14 to 0.24 depending on the orienta¬ 
tion of the cyclotron orbit m, in agreement with our 
result. 


V. FAR-INFRARED REFLECTIVITY 

The reflectivity for a sample thickness of d = 102 ^m 
is shown in the left panel of Fig. |S4[ Infrared-active 
phonon modes dominate the shape of R{uj) below about 
150 cm“^ (see also Refs. [S3l [S8l [Sol ISTol [Sm [ST^ [ST3l 
IS141IS151IS161 IS17| ). In the transparent frequency range 
above the phonons, the measured reflectivity shows in¬ 
terference fringes. For these measurements with polariza¬ 
tion parallel to the cleavage plane, we expect two phonon 
modes with symmetry |S17| . These modes were ob¬ 
served at 48 cm“^ and 98 cm“^ in Bi 2 Te 3 , at 61 cm“^ and 
134cm“^ in Bi 2 Se 3 |S17| . and at 62cm“^ and 120cm“^ 
in BiSbTeSe 2 [S8] . Since the spectral weight of the lower 
mode is much larger, the reflectivity data typically show 
one strong Reststrahlenband above about 50cm“^ with a 
weaker feature on its high-frequency side. In BiSbTeSe 2 , 
we find a strong Reststrahlenband with a weak shoul¬ 
der at about 124 cm“^, in reasonable agreement with the 
data reported by Reijnders et al. |S^ and by Post et al. 
[S8| . The strong Reststrahlenband is located close to the 
lower limit of our frequency range, which hinders the de¬ 
termination of the precise eigenfrequency. Compared to 
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FIG. S4: Reflectivity and dielectric function e(a;) in the far- 
infrared range. Left: Reflectivity for a sample thickness 
of d = 102 (solid lines) compared to Drude-Lorentz fits 
(dashed lines). Right: Dashed lines depict the dielectric func¬ 
tion e{uj) determined from the Drude-Lorentz fits of the reflec¬ 
tivity shown in the left panel. The Reststrahlenband carries 
a huge oscillator strength. Solid lines show e{Lj) as obtained 
from the interference fringes in the transparent range. 


Bi 2 Te 2 Se, the features are much broader in BiSbTeSe 2 , 
which can be attributed to the Bi/Sb and Te/Se disorder. 


VI. DIELECTRIC FUNCTION e{uj) 


For w —0, e{uj) is dominated by the huge oscillator 
strength of the strong phonon mode, see right panel of 
Fig. S4 A Drude-Lorentz fit yields e(a; —)■ 0) of the 


order of 200. A more precise determination of this value 
is difficult because the relevant phonon is close to the 
lower frequency limit of our data, see left panel of Fig.|S4| 
Similar values can be derived from results reported in the 
literature. In Bi 2 Te 2 Se, an analysis of the infrared-active 
phonon modes yields one strong phonon mode with a 
contribution of about 230, six more modes with a total 
contribution of about 10, and a value of £« 35-40 above 
the phonons [SlOj . In total, this yields s{uj — >■ 0)«275- 
280 in Bi 2 Te 2 Se. In Bi 2 Se 3 , a Drude-Lorentz fit yields a 
phonon contribution of about 170 and a value of £ « 30 
above the phonons, adding up to £(a; —>■ 0)^200 isnj. 


VII. FINITE-SIZE EFFECTS OF THE 
NUMERICAL DATA 


At T = 0 and for perfect compensation, AT =1, screen¬ 
ing is weakest and one expects the largest finite-size ef¬ 
fects in the numerical simulations. In this regime the typ¬ 
ical distance of puddles scales with (A/Ac)^. In Fig. S5 
we show ppA/Ec, as function of the inverse linear sys¬ 
tem size, where Pp denotes the fraction of dopants con¬ 
tributing to electron-hole puddles (see Fig. 5 of main 
text). While small system-size effects are clearly visi- 
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ble, all finite-size errors are well below 5%. Note that 
we expect a universal value of Pp IS./Ec for IS /> oo. 
For the shown values of A/Fic > 9 this limit is reached 
with only small corrections: PpA/Ec slightly grows for 
smaller values of A/Ec- 



FIG. S5: To illustrate finite-size effects of our numerical re¬ 
sults, ppA/Ec (see Fig. 7 of the main text) is plotted as func¬ 
tion of the inverse of the linear size L of the simulated box 
for perfect compensation K = 1 and T — 0 for A/Ec = 9, 12, 
and 15. The number of dopants is 2 x L^, reaching more 
than 400.000 for L = 60. Error bars represent a single stan¬ 
dard deviation of the mean value obtained by averaging over 
200 — 500 impurity configurations. All finite size effects are 
well below 5%. 
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